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We present a first-principles study of the relationship between stress, temperature and electronic 
properties in piezoelectric ZnO. Our method is a plane wave pseudopotential implementation of 
density functional theory and density functional linear response within the local density approxi- 
mation. We observe marked changes in the piezoelectric and dielectric constants when the material 
is distorted. This stress dependence is the result of strong, bond length dependent, hybridization 
between the O 2p and Zn 3d electrons. Our results indicate that fine tuning of the piezoelectric 
properties for specific device applications can be achieved by control of the ZnO lattice constant, 
for example by epitaxial growth on an appropriate substrate. 
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I. INTRODUCTION 

Zinc oxide (ZnO) is a tetrahedrally coordinated wide band gap semiconductor that crystallizes in the wurtzitc 
structure (Figure QJ). The lack of center of symmetry, combined with a large electromechanical coupling, result in 
strong piezoelectric properties, and the consequent use of ZnO in mechanical actuators and piezoelectric sensors. In 
addition, ZnO is transparent to visible light and can be made highly conductive by doping. This leads to applications in 
surface acoustic wave devices and transparent conducting electrodes. However the piezoelectric properties can change 
the characteristics of potential energy barriers to mobile charges at interfaces, and hence affect the carrier transport 
properties. The resulting piezoresistance is at times desirable, for example in ZnO-based metal-oxide varistors which 
can dissipate large amounts of power in short response times and are commonly found as electrical surge protectors 
M. However the detailed effects of piezoelectrically induced changes on the electrical behavior of ZnO have not yet 
been well characterized, and as ZnO finds increased application in electronic devices these effects will have large 
technological impact. 

In this paper we present a first-principles study of the strain dependence of the electrical properties of ZnO, using 
a plane wave pseudopotential (PWPP) implementation of density functional theory (DFT) within the local density 
approximation (LDA). We calculate and analyze the structural dependence of total energies, band structures and 
piezoelectric and dielectric constants. The principal result of our analysis is that the piezoelectric constants of ZnO 
are strongly dependent on conditions of stress and temperature, whereas the dielectric coefficients vary less strongly. 
Of particular interest is a comparison of the properties of ZnO, in which the Zn 3d bands are filled with those of 
, PbTi03, which has empty Ti 3d states. Both materials have anomalously large piezoelectric coefficients, but in 
PbTiOa the dielectric coefficients are also anomalously large, whereas ZnO behaves as an ordinary dielectric. Our 
hypothesis is that the strong O 2p - Zn 3d hybridization which has been previously noted in PbTi03 |^|| also occurs 
O ■ in ZnO, leading in both materials to strong strain-phonon coupling and consequently large piezoelectric coefficients. 
Thus a filled 3d band does not preclude a large peizoelectric response. The large dielectric response in PbTiOa has a 
similar origin. In ZnO however, the filled O 2p and Zn 3c? bands do not allow a large orbital response to an applied 
electric field, resulting in a normal dielectric response. 

The remainder of this paper is organized as follows. In Section [n] we summarize the results of earlier theoretical 



studies of ZnO. In Section III we describe the theoretical and computational methods used in this work. In Section IV, 
we describe our results for the static and response properties of bulk ZnO at its equilibrium lattice constant. In Section 
^ we investigate the dependence of the electronic structure and response functions on external stress or changes in 
temperature. Finally, in Section [vj we present our conclusions and discuss implications for materials growth and 
device design. 
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II. PREVIOUS THEORETICAL WORK 



First-principles studies of ZnO are computationally challenging. First, the wurtzite structure contains twice as many 
atoms per unit cell as the zincblende semiconductor structure. In addition, both oxygen and zinc are problematic 
atoms for the construction of pseudopotentials. In both cases the relevant valence electrons (O 2p and Zn 3d) 
have no lower-lying electrons of the same angular momentum to provide an effective repulsive potential from the 
orthogonalization requirement p|. As a result they are tightly bound and require a large number of plane waves in 
their expansion. 

The first published band structure of ZnO (|] used the Green's function KKR method. This was followed by 
empirical pseudopotential calculations || in which the Zn 3d electrons were placed in the core. In addition 
to preventing assessment of the Zn 3d contribution to bonding properties, this approach was later shown to give 
unsatisfactory results || . 

Schroer and co-workers Q circumvented the problem of a large plane wave basis set by combining the use of 
pseudopotentials with localized Gaussian basis sets containing orbitals of s, p, d and s* symmetry. Using this basis 
they compared the results of LDA calculations using Zn 2+ (in which the 3d electrons are contained in the core) and 
Zn 12+ (in which the 3d electrons are valence electrons) pseudopotentials, and found that the Zn 12+ pseudopotential 
gave results in good agreement with experiment. For example, their LDA energy minimum volume was 0.6 % below the 
experimental value and their calculated band structure was in reasonable agreement with angle resolved photoemission 
measurements. There was a slight discrepancy in the position of the d bands which they attributed to inadequacy of 
the local density approximation in describing these strongly correlated bands. 

Dal Corso et al. || avoided the use of pseudopotentials entirely by using the all electron full potential linear 
augmented plane wave (FLAPW) method. They calculated the piezoelectric and polarization properties of ZnO 
within the LDA, with the purpose of determining the origin of the unusually strong piezoelectric response in ZnO. 
The principal result of their work was that the contribution to the macroscopic polarization tensor from the relative 
displacement of the sublattices was large, and only partly canceled by the electronic "clamped-ion" contribution, 
leading to a large net piezoelectric polarization. (In contrast, in zincblende semiconductors these terms are of similar 
magnitude and opposite sign, resulting in a small piezoelectric polarization.) As in Ref. S, their LDA volume slightly 
underestimated the experimental value, and their Zn d bands were around 4 eV higher in energy than observed in 
photoemission. 

Hartree-Fock calculations have been used successfully to determine the stability of, and transitions between, different 
phases of ZnO |l(J and to calculate the (1010) surface reconstruction for the wurtzite phase p3]| . The Hartree-Fock 
approximation also gives an incorrect energy position for the Zn 3d bands - this time around 2 eV too low. 



III. COMPUTATIONAL TECHNIQUES 
A. Pseudopotential construction 

The calculations described in this work were performed using a plane wave pseudopotential implementation |L^ | 
of density functional theory [ p^[ within the local density approximation. Plane wave basis sets offer many advan- 
tages in total energy calculations for solids, including completeness, an unbiased representation, and arbitrarily good 
convergence accuracy. They also allow for straightforward mathematical formulation and implementation, which is 
invaluable in the calculation of Hellmann-Feynman forces |l4|| and in the density functional theory linear response 
calculations employed here p5[ . 

However plane wave basis sets necessitate the use of pseudopotentials to model the electron-ion interaction, in 
order to avoid rapid oscillations of the valence wavefunctions in the region around the ion cores. The difficulties 
associated with applying the pseudopotential method to tightly bound d-electrons, which might be expected to require 
a prohibitively large number of plane waves to expand their pseudopotentials, were mentioned above. An earlier study 
of cubic ZnS Jl6{ using the smooth Trouiller-Martins pseudopotentials, required plane waves up to 121 Ry in energy to 
achieve an energy convergence of 0.05 eV. Although feasible for a bulk calculation for the zincblende structure (with 
only two atoms per unit cell) such a large energy cutoff is undesirable for larger unit cells, such as that of the wurtzite 
structure, or those required for calculation of surface properties. In this work, we use the optimized pseudopotentials 
developed by Rappe et al. jl7|], which allow us to reduce the required energy cutoff to 64 Ry without compromising 
accuracy or transferability. Optimized pseudopotentials minimize the kinetic energy in the high Fourier components 
of the pseudo wavefunction, leading to a corresponding reduction in the contribution of high Fourier components in 
the solid. 



2 



For both Zn and O we constructed non-relativistic optimized pseudopotentials. The oxygen pseudopotentials were 
generated from a 2s 2 2p 4 reference configuration with core radii, r C: of f .5 a.u. for both s and p orbitals. They were 
then optimized using 4 and 3 basis functions with cutoff wave vectors, q c , of 7.0 and 6.5 a.u. for s and p orbitals 
respectively. q c determines the convergence of the kinetic energy with respect to the plane wave cutoff energy in 
reciprocal space calculations. These oxygen pseudopotentials were used in earlier calculations for perovskite oxides 
jisf | fl9f and gave accurate results. The zinc pseudopotentials were constructed for a neutral Zn atom with reference 
configuration 3c£ 10 4,s 175 4p - 25 . r c values of 2.0, 1.4 and 1.4 a.u.s were used for d, s and p orbitals respectively, with q c 
values of 8.0, 7.0 and 8.0 Ry (giving a cutoff energy of 64 Ry.) The transferability of the pseudopotential was tested 
for a variety of +1 and +2 free Zn ions. The pseudo total energies and eigenvalues were in agreement with the all 
electron values to within 0.001 a.u.s. There was no improvement in transferability, or in agreement with all electron 
calculations for bulk systems, on inclusion of non linear core corrections p0| . All pseudopotentials were put into 
separable form pl[ using one projector for each angular momentum. For both Zn and O the I = 1 component was 
chosen as the local potential. The absence of ghost states was confirmed using the ghost theorem of Gonze, Kackell 
and Schemer ||. 

B. Density functional theory linear response 

We use density functional theory linear response (DFT-LR) to obtain the quadratic couplings between homogeneous 
strain, internal displacements of atoms and macroscopic electric field ]2^| . We use a variational formulation of DFT- 
LR |24[] in which the second derivative of total energy is minimized with respect to the first derivatives of Kohn-Sham 
wave functions with appropriate orthogonality constraints. This method avoids using any finite-difference formulae 
and yields dielectric or piezoelectric constants with a minimal number of calculations. 

To obtain the piezoelectric and dielectric constants, we calculate the DFT-LR of our system to two types of 
perturbations: (a) phonon (or atomic displacements) and (b) electric field. Using the first-order response wavefunction 



resulting from (a) in the Hellman-Feynman force formula 14 and in the stress formula P5j, we obtain the dynamical 



matrix and the coupling between phonons and strain respectively. Similarly, the response wavefunctions resulting 
from (b) are used to obtain the Born effective charges and the clamped-ion piezoelectric constants respectively. 

The symmetry of the wurtzite structure allows three independent piezoelectric constants (733,713,714) and two 
dielectric constants (£33,613). In the present work, we focus on 733,713 and £33. The piezoelectric constants 733 and 
713 give the polarization along the c-axis induced by strains 633 and en respectively. Equivalently, 733,713 give the 
stresses 033,011 induced on the unit cell by an electric field along the c-axis. The former relationship underlies the 
earlier work of Dal Corso et al || based on finite-difference formulae and geometric phase, and the latter is used in 
the present work based on DFT linear response. 

To obtain the piezoelectric and dielectric constants, we perform two DFT-LR calculations, one with phonon per- 
turbations corresponding to Zn and O displacements along the c-axis, and another with a perturbing electric field 
along the c-axis. The phonon perturbation allows us to obtain both the frequencies of the T-point phonons with 
z-polarization (which are proportional to the square roots of the respective force constants, y/K a ), and the strain- 
phonon couplings, L azz and L axx . The LR calculations with the electric field perturbation yield the optical dielectric 
constant £33, the Born effective charges Z a , zz , and the clamped-ion piezoelectric couplings 733 and 7J3. The dielectric 
and piezoelectric constants are then given by: 

,00 1 j„ Z a zz Z a zz 
e 33 = £33 + 47r 2^ ' 



733 = 733 + ^ 



^OL^ZZ^OL^Z 



^a^zz-^a^xx 



713 = 7l3 + 2^ ' 
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C. Technical Details 



All calculations were performed on Silicon Graphics 02 and Origin 200 systems using the conjugate gradient 
program CASTEP 2.1 (2^J|^] and our own related density functional linear response program |28|| . We used a plane 
wave cut off of 64 Ry, which corresponds to around 3000 plane waves per wavefunction in a a single ZnO wurtzitc 
unit cell. A 3x3x2 Monkhorst-Pack |29| grid was used, leading to six k-points in the irreducible Brillouin Zone for 
the high symmetry structures, and a correspondingly higher number for distorted structures with lower symmetry. 
The exchange correlation was parameterized using the Perdew-Zunger parameterization |30| of the Ceperley-Alder 
potential . 



IV. UNSTRAINED ZNO; NEW RESULTS AND COMPARISON WITH LAPW CALCULATIONS 

Before presenting our results for strained ZnO, we first discuss bulk, unstrained ZnO, and compare the results 
of our plane wave pseudopotential calculations with published theoretical and experimental data. A good test of 
the transferability of the pseudopotentials is that they predict the same minimum energy structure as all-electron 
calculations which use the same approximation for the exchange-correlation functional. The minimum energy volume 
obtained by the earlier all-electron LDA calculation (Ref. ||) is 45.89 A 3 . Ref. || also showed that the LDA - ratio 
and u value are the same as the corresponding experimental quantities to within 0.5 mRy, and that the energy surface is 
rather flat around the minimum in the (it, ^) plane. The minimum energy volume obtained using our pseudopotentials 
is 47.6 A 3 (at the experimental £ ratio and u value), and that obtained by an earlier pseudopotential calculation using 
a Gaussian basis is 46.80 A 3 . The experimental volume is 47.90 A 3 and we believe the results obtained using our 
method are within LDA errors. 



A. Band structure analysis 

In Figure ^| we show the calculated band structure for ZnO at the experimental unit cell volume, (a = 3.2595 A, 
c = 5.2070 A and u = 0.3820), which is approximately equal to the LDA minimum volume for our pseudopotentials. 
The top of the valence band is set to eV. The band structure is indistinguishable from that presented in Ref. |jj 
which was calculated using the FLAPW technique. 

The narrow bands between around -6 eV and -5 eV derive largely from the Zn 3d orbitals and are completely 
filled. The broad bands between around -5 eV and eV are from the O 2p orbitals and are again completely 
filled. The Zn 4s band is broad and unoccupied, ranging between 1 and 7 eV above the band gap. Figure shows 
the same band structure along T to A with the symmetry labels within the Cq v group added. 1 is the totally 
symmetric representation, 2 is antisymmetric with respect to Cg and Ci rotation and o~ v reflection, 5 and 6 are the 
doubly degenerate representations, with 5 being antisymmetric with respect to C§ and C3 rotation, and 6 being 
antisymmetric with respect to C3 and C2 rotation. We see that interactions between O 2p and Zn 3c? orbitals are 
allowed by symmetry at T and A, and along the adjoining A line. The Zn s orbitals can interact with O 2p and Zn 
3d bands of 1 and 3 symmetry. 

In order to quantify the interactions between the various orbitals we perform a tight-binding analysis along the T 
to A direction of the Brillouin zone. Tight-binding parameters are obtained by non-linear-least-squares fitting ]3^ ] to 
the calculated ab initio energies at the high symmetry T and A points, and at 19 points along the A axis. A good 
tight-binding fit (rms deviation = 0.11) is obtained when only nearest neighbor interactions between O 2p and Zn 3d 
bands, and O 2p and Zn 4s bands are included in the fit. Additional small Zn 3d - Zn 3c? interactions are needed to 
produce dispersion in the upper e g Zn 3d band along this symmetry axis. The tight-binding parameters which we 
obtain are given in Table || (in the column labeled 'structure 1'), and the tight-binding band structure is compared 
with the ab initio band structure in Figure ^. 

The two sets of Zn-0 parameters correspond to the two Zn-0 distances, the shorter ( r\ = 1.974 A ) being the 
separation of the Zn and O atoms lying directly above each other along the c direction, and the longer ( r% = 1.989 
A ) joining Zn and O atoms in adjacent c-oriented 'chains'. The Zn 3d - Zn 3c? interactions are small. The Zn 4s 
- O 2p interactions are large, and show the expected increase when the Zn-0 spacing is decreased. The Zn 3 c? - 
O 2p interactions are also large, and are quite different (even changing sign) for the two different Zn - O atom pair 
types. In fact both the a and n components are larger for the in-plane pairs which have the larger Zn-0 spacing. 
This distance dependence is unusual for tight-binding parameters, and suggests that the nature of the Zn 3c? - O 2p 
hybridization is different for the Zn-0 pairs lying in the c axis chains, than for the basal plane Zn-0 pairs. 
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B. Piezoelectric and dielectric properties 



In Table E, we present our results for piezoelectric and dielectric compliances and a comparison with those obtained 
in previous FLAPW calculations |9) and experiments Q . The agreement between the computational results is good, 
but the scatter in experimental results indicates that the piezoelectric constants may be sensitive to experimental 
conditions. We will investigate this issue in the next section. We also report in Table | the frequencies of TO phonons 
with c polarization at the T point, and note that only one of the TO phonons (395 cm -1 ) is IR-active; this is the only 
phonon which contributes to piezoelectric and dielectric constants 733,713 and £33, studied in this work. Our results 
also confirm the conclusion of Ref. II that the piezoelectric constants of ZnO are dominantly contributed by phonons 
(internal strain). The dielectric constant, on the other hand, has roughly equal contributions from electrons (e°°) and 
phonons. 

V. EFFECTS OF STRAIN AND TEMPERATURE ON THE ELECTRONIC PROPERTIES OF BULK ZNO 



Next we investigate the effects of strain on the static electronic properties of bulk ZnO. We compare two different 
strained structures with the equilibrium structure at the experimental lattice constant. First we simulate application 
of a homogeneous in-plane compression by reducing the a lattice constant by 2%, while increasing the c lattice 
constant correspondingly to maintain the same total volume. The value of u is held at 0.3820 of the c axis. Second 
we investigate the effect of changing the Zn-0 separation along the c axis, u, by increasing it by 5% to 0.4011, while 
retaining the equilibrium a and c values. 



Application of the in-plane compression to reduce the a lattice constant by 2% increases the total energy by 0.1 
eV, and creates Hellman-Feynman forces on the atoms in the z direction of ±0.59t. The band structure of the 

compressed structure is shown in Figure pi We observe a broadening of the bands, as expected from the increased 
overlap between the orbitals in the basal plane. Note in particular the broadening in the Zn 3d bands. In spite of the 
fact that the d bands are narrow, they are very sensitive to strain and bond length as a result of p — d hybridization 
processes. 

In addition to an overall band broadening, there are a number of significant details in the band structure for this 
compressed structure. First, the O 2p bands with symmetries 1 and 3 shift down relative to the uppermost O 2p 
bands at the T point. This leads to a reduction in the density of states at the Fermi level, and to an overlapping of 
the low energy part of these bands with the Zn 3d bands. The ordering of the lower Zn 3d bands is reversed compared 
with those in the unstrained structure, and the Zn 4s bands shift slightly down in energy, resulting in a smaller band 
gap. 

These observations are consistent with a tight-binding fit along the r to A line, the results of which are given in 
Table [jl] (structure 2). Again, a good (root mean square deviation = 0.11) tight-binding fit is obtained using the 
limited interaction set described above. Again the O 2p - Zn 4s hybridization is large and shows the expected variation 
with bond length. The dependence of the (also large) p — d parameters on distance does not follow a straightforward 
pattern. Note that the anomalously low value for the Zn 4s energy is the result of an attempt by the fitting package 
to reproduce the down shift in the O 2p band of 1 symmetry (which also contains a significant Zn 4s component) 
within our limited basis set. This energy value shifts up to a more physical positive number if additional interactions 
are included in the basis. 



Change in the u value by 5 % increases the total energy by only 0.07 eV compared with the equilibrium structure, 
and introduces Hellman Feynman forces of ±0.59 eV in the z direction on the ions. The new band structure is shown 
in Figure |6|, and is qualitatively very similar to that of the undistorted structure. 

A tight-binding fit using the interaction set described above produces an rms deviation of 0.13, and the parameters 
listed in Table |l|. Again the O 2p - Zn 4s interaction is largest for the smallest bond length as expected, but the p — d 




A. Band structures 



1. In-plane compression 



2. Change in u value 
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hybridizations have a more complicated distance dependence. In this case, the O 2p bands with 1 and 3 symmetries 
are shifted down a small amount at T, but not as much as in the previous structure. As a result they do not overlap 
with the Zn 3c? bands. The Zn 4s bands are not shifted down relative to their position in the unstrained case. 



B. Piezoelectric properties 

1. Stress dependence 

Finally we explore the piezoelectric and dielectric properties of strained ZnO in detail. As a result of anharmonic 
couplings between different phonons and strain, the properties of ZnO are strongly structure sensitive. There is a 
large region in the phase space of structural configurations which is low in energy and therefore contributes to the 
properties of ZnO. To obtain energies, polarization and other properties for each of these configurations from first- 
principles is inefficient and impractical. We choose instead to use a model energy functional that captures the physics 
of the low-energy structural excitations of ZnO. 

To keep the model simple, we restrict our analysis to the subspace of degrees of freedom that preserve the symmetry 
of the wurtzite structure { x = e xx = e yy , z = e zz , u }. Justification for this simplification rests on the assumption 
(verifiable through ab initio calculations) that most of the low energy configurations preserve wurtzite symmetry and 
most of the symmetry breaking distortions are described well at harmonic order and can be integrated out. This is in 
contrast with earlier work on structural phase transitions |Q where the symmetry-breaking distortions were retained 
in the subspace. 

We write the model energy functional as a symmetry-invariant Taylor expansion in atomic displacements (or normal 
mode degrees of freedom u) and strains (x, z) defined above. Including the lowest order coupling with external stress 
a a and electric field along c-axis E, 

E tot {x, z, u, a a ,E) = ^Ku 2 + Au 3 + Bu A + ^[C x x 2 + C z z 2 + C xz xz] 
+D x x 3 + E x x A + D z z 3 + E z z A + V F a r a u 2 

a 

+ ^ G a/3r a r l3 u + ^ Ha/3 

a,0 a, (3 

- fl(2a x x + a z z + 2j^ 3 xE + j° 3 zE) - ZuE - —e^E 2 , (1) 

An 

where the 7°'s and are clamped (or electronic) piezoelectric and dielectric constants, and r a is strain x or z. K, 
A ... H, 7 and e are the harmonic and anharmonic coupling parameters including force, elastic and mixed coupling 
constants. These parameters have been determined from DFT total energy and linear response calculations. 

The equilibrium state of ZnO under the application of external stress or electric field at zero kelvin is obtained by 
minimizing the total energy E to t with respect to the structural parameters x, z and u: 

dE to t _ q dEtot _ q dEtot _ q 
dx ' dz ' du 

For the equilibrium structure, the static dielectric and piezoelectric constants are calculated using the expressions in 



Section III (also in Ref. |23J), with various coupling parameters being dependent on the structure. The spontaneous 
polarization is calculated using the expression 

p _ dE tot 



dE 

ZnO is grown in the form of thin solid films on sapphire and there are always interfacial stresses in the films. We 
use our model to study properties of ZnO as a function of stress in the basal ab plane. In Fig. ^, we show our results 
for the structural parameter u, dielectric and piezoelectric response of ZnO as the applied stress a — a xx — o~ yy is 
varied from -1 to 1 GPa. 

We find that parameter u changes by about 1 percent in this range of basal stress. The dielectric constant £33 
monotonically increases by about 2 to 3 percent. In contrast, the piezoelectric constants are rather sensitive to stress 
changing by about 15 to 30 %. With access to different contributions to these constants in our model, we discover 
that most of the dependence of these compliances on stress is due to the phonon contribution. 
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2. Temperature Dependence 



ZnO is also a pyroelectric material, characterized by a temperature dependence of the polarization. This has 
technological relevance because it leads to the widespread use of ZnO in infra-red detectors. To explore various 
stress and electric field-dependent properties of ZnO at finite temperature, we obtain free energy functional using a 
local harmonic model |?4j for entropy. In this model, entropy is calculated treating phonons harmonically for given 
structural parameters. In the present work, we include the optical phonons that are polarized along z-axis. With 
these approximations, the free energy is: 

G(x, z, u, a x ,a z ,E, T) = E tot (x, z, u, a x , a Z7 E) + 3k B TLog( — — ) (3) 

kbT 

where the K^s are the harmonic force constants of the three T— phonons with z-polarization for given values of 
structural parameters. Due to the anharmonic terms included in the energy expansion Eqn (1), Ki's (hence the 
T-dependent part of the free energy) depend on the structure. The equilibrium state of ZnO under the application 
of external stress or electric field at finite temperature is then obtained by minimizing the free energy G with respect 
to the structural parameters x, z and u: 

— = — = 0— = (4) 

dx ' dz ' du 



Again, the dielectric and piezoelectric constants are calculated using the expressions in Section |I|, and the sponta- 
neous polarization is now calculated using the expression 

dG 
BE' 

We investigate the dependence of various properties of ZnO on temperature. In Fig. |t], we display results for 
structural parameters, dielectric and piezoelectric constants, and polarization for a range of temperatures from to 
450 K. We find that the structural parameters such as the bond-length u and volume change only by about 0.3 % 
from zero kelvin to room temperature. 

The pyroelectric constant, which we obtain from our calculated results for polarization, is 20 /xC/m 2 /K, compared 
with the experimental value of 9.4 \xC jm 2 j K f33|| . Considering the simplicity in our treatment of temperature through 
the model entropy function, we find this agreement quite encouraging. In particular, we point out that, while the 
nonpolar TO phonons with z-polarization have been omitted from the expansion of energy, they have been included 
in the entropy term, where we use the expression of entropy J34| ] treating phonons harmonically, consistent with our 
energy expansion. 

To estimate the effect of these non-polar phonons on the pyroelectric constant, we omitted their contribution to 
entropy and found a pyroelectric constant of about 40 nC/m 2 /K (almost doubled). If we assume that phonons 
should generally suppress the pyroelectric constant, the discrepancy between theory and experiment is likely due to 
our omission of phonons with x and y-polarization from the entropy expression. 

The linear dielectric response 633 of ZnO is quite sensitive to temperature, changing by about 4% in the temperature 
range considered. The piezoelectric response, on the other hand, is very sensitive to temperature changing by about 
20%. This should be an important consideration in designing piezoelectric devices for operation at room temperature. 



C. Discussion 



It is clear from our results that the piezoelectric response of ZnO is strongly sensitive to both temperature and 
stress, changing by up to 30 % over the range of parameters considered. This dependence arises from the changes 



in structural parameters (manifested through the phonon contribution). We saw in Section III that the phonon 
contribution to the piezoelectric constants arises from the coupling of phonons with strain, L, and Born effective 
charge, Z, 

LZ 

^phonon — 7~7 , 
K 

K being the force constant. In Fig. ^, we show how K and L for the polar TO phonon change with temperature. 
While K changes by only 10 %, the coupling with strain, L, changes by about 25 %. The Born effective charge Z 
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(which describes the coupling of this phonon with electric field) is not found to vary much with temperature. The 
large temperature dependence of the piezoelectric response arises predominantly from that of the coupling of phonon 
with strain and its force constant. Since only the latter contributes to the dielectric constants, dielectric properties 
are less sensitive to s tructure, stress or temperature. 

In the Section V A , we found the hybridization between Zn d orbitals and O p orbitals to be sensitive to structural 
parameters. The same hybridization was found to be the cause of anomalous Born effective charges in ferroelectric 
materials such as PbTiC>3 ||. While the d orbitals of the transition metals in perovskite ferroelectrics are formally 
unoccupied, those in Zn are fully occupied, leading to normal effective charges. The coupling of phonons with strain, 
however, is large and structure dependent irrespective of the occupancy of d orbitals. 



VI. SUMMARY 



In summary, we have calculated the electronic and atomic structure of ZnO from first-principles, and analyzed the 
nature of the bonding using the tight-binding method. Wc find that hybridization between the Zn d orbitals and 
O p orbitals is strongly structure dependent. Using DFT linear response, we have obtained the phonon frequencies, 
dielectric and piezoelectric constants of ZnO at zero temperature, and have shown that phonons (internal strain) 
have the dominant contribution to piezoelectricity in ZnO. ^From DFT linear response and total energy calculations, 
and a simple model for vibrational entropy, we have constructed an ab initio free energy functional for ZnO to study 
its properties at finite temperature and under applied stress. Our results show that the piezoelectric properties of 
ZnO are strongly dependent on both temperature and stress. This clearly has implications for the design of devices 
intended to operate at room temperature, or under stressed conditions. By analyzing various physical contributions, 
we have found that this is primarily due to the coupling between phonons and strain. The O 2p - Zn 3d hybridization 
is the cause of the large magnitude and sensitivity of this coupling. 
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property 


Ref. [3] 


Present Work 


Experiment 33j 


73°3 (C/m 2 ) 
733 (C/m 2 ) 
7?3 (C/m 2 ) 
7i3 (C/m 2 ) 

e 33 
£33 

TO phonon (cm" 1 ) 


2.05 
-0.58 
1.21 
0.37 
-0.51 


2.07 
-0.73 
1.30 
0.31 
-0.66 
4.39 
8.75 

544.919, 395.349, 258.519 


2.10 
1.0-1.5 
-0.36 to -0.62 



TABLE I. Comparison between results of FLAPW calculations (Ref. |9|), pseudopotential calculations (this work) and 
experimental results (Ref. []33J) for piezoelectric constants and related properties of ZnO. 



parameter 


structure 1 


structure 2 


structure 3 


Eo2p 


-1.262 


-1.062 


-1.308 


Eznis 


1.235 


-0.274 


1.409 


Eznid 


-5.144 


-5.396 


-5.159 


Vo2p-Znis 


2.584 


2.936 


2.823 


V 2 

V 02p-Zn4.s 


2.386 


2.517 


2.279 


V (02p-Zn3d)cr 
V 2 

v (02p-Zn3d)cr 
V 1 

v (02p-Zn3d)Tv 
V 2 

v (02p-Zn3d)Tv 
V(Zn3d-Zn3d)<T 


0.625 
-0.856 

1.292 
-1.711 

0.169 


0.829 
-1.343 

1.444 
-1.365 

0.124 


0.591 
-0.761 

1.294 
-1.802 

0.170 


V(Zn3d-Zn3d)TT 


-0.006 


0.001 


-0.014 


V(Zn3d-Zn3d)S 


0.053 


0.034 


0.060 



TABLE II. Tight-binding parameters (in eV) for ZnO obtained by non-linear-least-squares fitting to the ab initio eigenvalues 
along r to A. E indicates an orbital energy, and V an inter-atomic transfer integral. The transfer integrals with the superscript 
'1' are between the closest nearest neighbor Zn-O pairs, and those with the superscript '2' are between the nearest neighbors 
with the larger separation. Only the parameters listed in the table were allowed to be non-zero in the fitting procedure. 
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FIG. 2. Calculated band structure of ZnO. 
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FIG. 4. Comparison of tight-binding and ab initio band structures in ZnO. 



13 




r A H K r 



FIG. 5. Calculated band structure of strained ZnO. 
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r A H K r 

FIG. 6. Calculated band structure of ZnO with the u value increased by 5 %. 
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FIG. 7. Dependence of the change in structural parameter An (a), volume AV/V (b) and spontaneous polarization AP 
(c), dielectric constant £33 (d), piezoelectric constants 733 (e) and 713 (f) of ZnO on the temperature. 
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FIG. 8. Dependence of the change in structural parameter Am (a), dielectric £33 (b), piezoelectric constants 733 (c) and 713 
(c) of ZnO on the applied stress in the basal plane a = a xx = o yy . 
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FIG. 9. Dependence of (a) the force constant and (b)coupling with strain of the polar TO phonon as a function of 
temperature. 
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